home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / newmat08.zip / TMTF.CPP < prev    next >
C/C++ Source or Header  |  1995-01-11  |  4KB  |  140 lines
  1.  
  2. //#define WANT_STREAM
  3. #define WANT_MATH
  4.  
  5. #include "include.h"
  6.  
  7. #include "newmatap.h"
  8.  
  9. void Print(const Matrix& X);
  10. void Print(const UpperTriangularMatrix& X);
  11. void Print(const DiagonalMatrix& X);
  12. void Print(const SymmetricMatrix& X);
  13. void Print(const LowerTriangularMatrix& X);
  14.  
  15. void Clean(Matrix&, Real);
  16.  
  17.  
  18. static void SlowFT(const ColumnVector& a, const ColumnVector&b,
  19.    ColumnVector& x, ColumnVector& y)
  20. {
  21.    int n = a.Nrows();
  22.    x.ReDimension(n); y.ReDimension(n);
  23.    Real f = 6.2831853071795864769/n;
  24.    for (int j=1; j<=n; j++)
  25.    {
  26.       Real sumx = 0.0; Real sumy = 0.0;
  27.       for (int k=1; k<=n; k++)
  28.       {
  29.      Real theta = - (j-1) * (k-1) * f;
  30.      Real c = cos(theta);  Real s = sin(theta);
  31.      sumx += c * a(k) - s * b(k);
  32.      sumy += s * a(k) + c * b(k);
  33.       }
  34.       x(j) = sumx; y(j) = sumy;
  35.    }
  36. }
  37.  
  38. static void test(int n)
  39. {
  40.    Tracer et("Test FFT");
  41.  
  42.    ColumnVector A(n), B(n), X, Y;
  43.    for (int i=0; i<n; i++)
  44.    {
  45.       Real f = 6.2831853071795864769*i/n;
  46.       A.element(i) = fabs(sin(7.0*f) + 0.5 * cos(19.0 * f)) + (Real)i/(Real)n;
  47.       B.element(i) = fabs(0.25 * cos(31.0 * f)) + (Real)i/(Real)n;
  48.    }
  49.    FFT(A, B, X, Y);
  50.    FFTI(X, Y, X, Y);
  51.    X = X - A; Y = Y - B;
  52.    Clean(X,0.000000001); Clean(Y,0.000000001); Print(X); Print(Y);
  53. }
  54.  
  55. static void test1(int n)
  56. {
  57.    Tracer et("Test RealFFT");
  58.  
  59.    ColumnVector A(n), B(n), X, Y;
  60.    for (int i=0; i<n; i++)
  61.    {
  62.       Real f = 6.2831853071795864769*i/n;
  63.       A.element(i) = fabs(sin(7.0*f) + 0.5 * cos(19.0 * f)) + (Real)i/(Real)n;
  64.    }
  65.    B = 0.0;
  66.    FFT(A, B, X, Y);
  67.    B.CleanUp();                 // release some space
  68.    int n2 = n/2+1;
  69.    ColumnVector X1,Y1,X2,Y2;
  70.    RealFFT(A, X1, Y1);
  71.    X2 = X1 - X.Rows(1,n2); Y2 = Y1 - Y.Rows(1,n2);
  72.    Clean(X2,0.000000001); Clean(Y2,0.000000001); Print(X2); Print(Y2);
  73.    X2.CleanUp(); Y2.CleanUp();  // release some more space
  74.    RealFFTI(X1,Y1,B);
  75.    B = A - B;
  76.    Clean(B,0.000000001); Print(B);
  77. }
  78.  
  79. static void test2(int n)
  80. {
  81.    Tracer et("cf FFT and slow FT");
  82.  
  83.    ColumnVector A(n), B(n), X, Y, X1, Y1;
  84.    for (int i=0; i<n; i++)
  85.    {
  86.       Real f = 6.2831853071795864769*i/n;
  87.       A.element(i) = fabs(sin(7.0*f) - 0.5 * cos(19.0 * f)) + (Real)i/(Real)n;
  88.       B.element(i) = fabs(0.25 * cos(31.0 * f)) - (Real)i/(Real)n;
  89.    }
  90.    FFT(A, B, X, Y);
  91.    SlowFT(A, B, X1, Y1);
  92.    X = X - X1; Y = Y - Y1;
  93.    Clean(X,0.000000001); Clean(Y,0.000000001); Print(X); Print(Y);
  94. }
  95.  
  96. void trymatf()
  97. {
  98.    Tracer et("Fifteenth test of Matrix package");
  99.    Exception::PrintTrace(TRUE);
  100. //   cout << "\nFifteenth test of Matrix package\n";
  101. //   cout << "\n";
  102.  
  103.    ColumnVector A(12), B(12);
  104.    for (int i = 1; i <=12; i++)
  105.    {
  106.       Real i1 = i - 1;
  107.       A(i) = .7
  108.        + .2 * cos(6.2831853071795864769 * 4.0 * i1 / 12)
  109.        + .3 * sin(6.2831853071795864769 * 3.0 * i1 / 12);
  110.       B(i) = .9
  111.        + .5 * sin(6.2831853071795864769 * 2.0 * i1 / 12)
  112.        + .4 * cos(6.2831853071795864769 * 1.0 * i1 / 12);
  113.    }
  114.    FFT(A, B, A, B);
  115.    A(1) -= 8.4; A(3) -= 3.0; A(5) -= 1.2; A(9) -= 1.2; A(11) += 3.0;
  116.    B(1) -= 10.8; B(2) -= 2.4; B(4) += 1.8; B(10) -= 1.8; B(12) -= 2.4;
  117.    Clean(A,0.000000001); Clean(B,0.000000001); Print(A); Print(B);
  118.  
  119.  
  120.  
  121.    test(2048);
  122.    test(2000);
  123.    test(27*81);
  124.    test(2310);
  125.    test(49*49);
  126.    test(13*13*13);
  127.    test(43*47);
  128.    test1(98);
  129.    test1(100);
  130.    test1(2048);
  131.    test1(2000);
  132.    test1(35*35*2);
  133.    test2(13);
  134.    test2(12);
  135.    test2(9);
  136.    test2(16);
  137.  
  138. //   cout << "\nEnd of Fifteenth test\n";
  139. }
  140.